#determine current working directory
getwd()
#read in local file
tadpoles <- read.table("C:/Users/jmweiss/Documents/ecol 563/tadpoles.csv", header=T, sep=',')
#this yields an error because the extension is missing in the file name
tadpoles <- read.table("C:/Users/jmweiss/Documents/ecol 563/tadpoles", header=T, sep=',')
#The read.csv function reads in comma-delimited files directly
tadpoles <- read.csv("C:/Users/jmweiss/Documents/ecol 563/tadpoles.csv")
#examine beginning and end of the file
head(tadpoles)
tail(tadpoles)
#look at first six observations
tadpoles[1:6,]
#look at first six observations and first two variables
tadpoles[1:6,1:2]
 
#You can read in files from web sites but not if they're https URLs
junk <- read.csv("https://sakai.unc.edu/access/content/group/3d1eb92e-7848-4f55-90c3-7c72a54e7e43/public/data/tadpoles.csv")
 
#To read in https files we need the RCurl package and a two-step process
library(RCurl)
#This works in Mac OSX
myfile<-getURL("https://sakai.unc.edu/access/content/group/3d1eb92e-7848-4f55-90c3-7c72a54e7e43/public/data/tadpoles.csv")
#This works in Windows--need two more arguments
myfile<-getURL("https://sakai.unc.edu/access/content/group/3d1eb92e-7848-4f55-90c3-7c72a54e7e43/public/data/tadpoles.csv", ssl.verifyhost=F, ssl.verifypeer=F)
#Both Mac OSX and Windows: Finally convert the file into an internal R format
tadpoles1 <- read.csv(textConnection(myfile))
 
#You can also read in files interactively
temp.dat <- read.csv(file.choose())
 
#You can read in files from local folders using a forward slash between folders
tadpoles <- read.csv("C:/Users/jmweiss/Documents/ecol 563/tadpoles.csv")
# or by using a double back-slash (Windows). Single back slashes are not allowed
tadpoles <- read.csv("C:\\Users\\jmweiss\\Documents\\ecol 563\\tadpoles.csv")
 
# To access a variable inside a data frame we must also reference the data frame as part of the name
# To display the variable called response
tadpoles$response
# Or specify the first column
tadpoles[,1]
# Or specify the column named response
tadpoles[,"response"]
 
#determining the nature of R variables with the class function
class(tadpoles$response)
#because its values are characters treatment was read in as a factor with associated dummy variables
class(tadpoles$treatment)
contrasts(tadpoles$treatment)
#fac1 is a factor with associated dummy variables
class(tadpoles$fac1)
contrasts(tadpoles$fac1)
#fac2 is a factor with associated dummy variables
contrasts(tadpoles$fac2)
 
#fac3 is numeric because its values are numbers
class(tadpoles$fac3)
tadpoles$fac3
 
#we can convert a numeric variable to a factor with the factor function
tadpoles$fac3.f <- factor(tadpoles$fac3)
class(tadpoles$fac3.f)
contrasts(tadpoles$fac3.f)
 
#I overwrite the numeric variable version with the factor version
tadpoles$fac3 <- factor(tadpoles$fac3)
contrasts(tadpoles$fac3)
 
#analyze experiment as a one-way ANOVA model: assign results to an object
out0 <- lm(response~treatment, data=tadpoles)
#placing the output object on the right side
lm(response~treatment, data=tadpoles) -> out0
out0
#summary table
summary(out0)
#ANOVA table
anova(out0)
 
#analyze as a 3-way ANOVA model
out1 <- lm(response~fac1+fac2+fac3+fac1:fac2+fac1:fac3+fac2:fac3+fac1:fac2:fac3, data=tadpoles)
#same number of parameters in both cases, just different parameterization
coef(out0)
coef(out1)
length(coef(out0))
length(coef(out1))
 
#short-cut notation for fitting a three-factor interaction model
out1 <- lm(response~fac1*fac2*fac3, data=tadpoles)
out1 <- lm(response~(fac1+fac2+fac3)^3, data=tadpoles)
 
#sequential ANOVA table: Type I tests
anova(out1)
 
#ANOVA table with Type II tests
library(car)
Anova(out1)
 
#based on the Type II tests we can drop three factor interaction and fac1:fac3
out2 <- lm(response~fac1+fac2+fac3+fac1:fac2+fac2:fac3, data=tadpoles)
 
#alternatively we can drop terms from the full model using the update function
out2a <- update(out1, .~.-fac1:fac2:fac3-fac1:fac3)
 
#Type II tests and Type I tests disagree about the fac1:fac2 interaction
anova(out2)
Anova(out2)
 
#graphical display of interactions
library(effects)
#obtain means from each effect
effect('fac2:fac3', out2)
#graph the fac2:fac3 interaction
plot(effect('fac2:fac3', out2), multiline=T)
#graph the fac1:fac2 interaction
plot(effect('fac1:fac2', out2), multiline=T)
 
#swap the roles of fac1 and fac2 in the plot
plot(effect('fac1:fac2',out2), multiline=T, x.var='fac2', z.var='fac1')
 
#without multiline=T we get profiles graphed separately
plot(effect('fac1:fac2',out2), x.var='fac2', z.var='fac1')
 
#plot the 3-factor interaction
plot(effect('fac1:fac2:fac3',out1), multiline=T)
 
#Check to see if variance is constant across treatments
boxplot(response~treatment, data=tadpoles)
boxplot(response~treatment, data=tadpoles, horizontal=T)
 
#calculate mean response by treatment
tapply(tadpoles$response, tadpoles$treatment, mean)
#calculate mean response by treatment: remove missing values first
tapply(tadpoles$response, tadpoles$treatment, mean, na.rm=T)
#save results
tapply(tadpoles$response, tadpoles$treatment, mean, na.rm=T) -> my.means
#add the means to the boxplots
points(my.means, 1:12, col=2, pch=8)